Doublon production rate in modulated optical lattices 



(N 

o 

(N 

>v 
c3 



Akiyuki Tokuno, 1 Eugene Demler, 2 and Thierry Giamarchi 1 

1 DPMC-MaNEP, University of Geneva, 24 Quai Ernest- Ansermet CH-1211 Geneva, Switzerland. 
2 Department of Physics, Harvard University, Cambridge, MA 02138. 
(Dated: January 27, 2013) 

We study theoretically lattice modulation experiments with ultracold fermions in optical lattices. 
We focus on the regime relevant to current experiments when interaction strength is larger than the 
bandwidth and temperature is higher than magnetic superexchange energy. We obtain analytical 
expressions for the rate of doublon production as a function of modulation frequency, filling fac- 
tor, and temperature. We use local density approximation to average over inhomogeneous density 
for atoms in a parabolic trap and find excellent agreement with experimentally measured values. 
Our results suggest that lattice modulation experiments can be used for thermometry of strongly 
interacting fermionic ensembles in optical lattices. 
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I. INTRODUCTION 

Cold atoms provide a new platform in which one can 
explore long standing open questions of strongly corre- 
lated systems in condensed matter physics [H, HI. In 
particular, two-component Fermi mixtures in an opti- 
cal lattice provide an ideal realization of the fermionic 
Hubbard model, where two species of fermions - cor- 
responding to the spin 1/2 - interact with an on-site 
repulsion. This model is relevant for understanding 
properties of electrons in several classes of novel quan- 
tum materials including oxides and high-T c superconduc- 
tors [H, IH . There are currently many efforts to probe the 
low-temperature physics of such a model, with experi- 
ments already demonstrating the Mott insulating behav- 
ior expected in this model [B|, Q . 

An important feature of strongly correlated ultra- 
cold atoms is that traditional probes used in solid state 
physics are often not readily available. One thus needs 
to understand how experimental techniques appropriate 
for atomic ensembles can provide information on many- 
body states. In this paper we focus on understanding 
lattice modulation experiments with fermions in opti- 
cal lattices. The technique of lattice modulation was 
originally introduced for bosonic systems and absorbed 
energy was measured as a function of modulation fre- 
quency 0|- Measuring energy absorption, however, can 
not be done accurately enough for strongly interacting 
fermions. Thus an extension of this technique was pro- 
posed, [1] and implemented @, in which the number of 
doubly occupied sites created by the lattice modulation 
was measured. Recent experiments successfully reached 
the regime of weak perturbations in which the number of 
doublons created scales quadratically with the modula- 
tion amplitude (Fermi's golden rule) and the modulation 
pulse duration |9(. 

While theoretical understanding of such experiments 
with bosons is now relatively complete fToj - fl~2( | . the case 
of fermions turned out to be more challenging. The main 
difficulty is the presence of excitations at very different 
energy scales: high-energy charge excitations, so-called 



doublons and holons, that have energies set by the on- 
site repulsion U and fermion hopping strength J, and 
magnetic excitations that have energies of the order of su- 
perexchange energy J 2 /U. Understanding the interplay 
of charge and spin degrees of freedom in the Hubbard 
model is a long standing problem in condensed matter 
physics [III, [HI • In the special case of half filling and fully 
disordered spin states, analysis of lattice modulation ex- 
periments has been performed previously [TH [H| . How- 
ever, such analysis is not sufficient for quantitative com- 
parison to experiments which are done with systems in 
a parabolic potential that have a large number of atoms 
outside of the incompressible Mott plateaus. 

Furthermore such real-time dynamics at finite tem- 
perature as seen in current experiments is problematic 
for numerics. Monte Carlo simulations [l7| suffer from 
the problems of analytic continuation. Density-matrix 
renormalization group approaches, which can deal with 
real time dynamics, are so far limited to one dimension. 
Therefore, a comparison of dynamical quantities to ex- 
periments is a non-trivial theoretical subject. 

In this paper we develop another analytical approach 
to calculate finite-temperature dynamics, that is, the 
doublon production rate (DPR), based on the slave- 
particle technique [TH, [H, [l!| . This approach is particu- 
larly adapted to the paramagnetic phase of the Hubbard 
model, and can be applied to any filling of the band and 
finite temperatures whose region is relevant to ongoing 
experiments. It provides a remarkable agreement to the 
experiments and allows for potential extensions. 

This paper is organized as follows. We define the sys- 
tem Hamiltonian, and introduce the slave-particle rep- 
resentation in Sec. [TTJ In Sec. IIII1 the spectral functions 
of the slave particles are evaluated, and in addition the 
spectral function of the original fermionic atom is also ob- 
tained. In Sec. lIVl we proceed with the estimation of the 
doublon production rate based on the obtained spectral 
functions in Sec. lIIIl and the analytic formula of the DPR 
spectrum is given. In Sec. |Vl the obtained analytic for- 
mula is extended to an inhomogeneous system in a trap 
by the local density approximation (LDA), and a com- 



2 



parison with the experiment in Ref. |9] is implemented. 
Finally, Sec. IVIl is devoted to the summary. 



II. SLAVE PARTICLE REPRESENTATION 

We consider the Hubbard model Hq = Hk + H at with 
H K = -J 4a C j^ ( la ) 

H at = -/i Ujy + U njfUji, (lb) 
i,<y j 

where cj CT and rij a — <^j a Cja are, respectively, a creation 
and number operator of a spin-cr fermions at a jth site. 
We use the following slave-particle representation: 
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.,=6^-6^], (2) 



where ft,] and are, respectively, creation oper- 
ators of a slave boson of spin-er state (Schwinger bo- 
son), holon and doublon at a site j. They satisfy the 
(anti-) commutation relations, [biajb^,] — 5i.j5a t a> and 

{hi,hj} = {di,dj} = Sij. The enlarged Hilbert space 
is projected onto the physical one by the following con- 
straint at every site: 
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Equation © allows us to rewrite the Hamiltonian (UJ as 
Hk = J E [ F Ji( h i h J - d\dj) + (Aljdjhi + H.c.)] , (4) 
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where Fji = J2<? b ] a bia and A\j = b^b^-b^b^ mean the 
hopping of slave bosons and the creation of spin singlet 
pair, respectively. The local potentials of a slave boson, 
holon, and doublon are, respectively, defined as = Xj, 
e!? = /i + Xj , and = U — /i + Xj. The constraint J3]) is 
implemented via the Lagrange multiplier Xj . 



III. DIAGRAMMATIC ANALYSIS 

A. Atomic limit 

We start with the atomic limit ( J/U = 0). Then the ki- 
netic term ^ which describes the scattering among slave 
particles vanishes. Since the atomic Hamiltonian ([5]) is 
quadratic, the atomic propagators at jth site are easily 
obtained as 



,(0) 
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(6) 



where u> n and v n are the Matsubara frequency for bosons 
and fermions, respectively. Note that the atomic limit 
propagators of slave bosons are independent of the spin. 
This means that the atomic limit exhibits spin-incoherent 
paramagnetism. Hereafter we set % = 1. 

Let us suppose the mean-field (MF) to be deter- 
mined by the atomic limit. Namely, the self-consistent 
equation for Xj corresponds to the statistical average of 
the constraint ([3]) in the atomic limit: 



2&(e bj ) + /(e hj ) + /(e dj ) = l, 



(7) 



where the prefactor 2 comes from the spin degrees of free- 
dom. /(e) and 6(e) are, respectively, the Fermi and Bose 
distribution functions. One can expect that if the effect of 
the kinetic energy Hk is small, that is, at relatively high 
temperature compared to the kinetic energy, the valid- 
ity of this treatment should be guaranteed. We thus use 
the MF assumption for the Lagrange multiplier: Xj X. 
Simultaneously the local potentials are also replaced by 
the homogeneous ones: e* — » e x where x = b, h, d. Cor- 
responding to the MF treatment for the Lagrange multi- 
plier, the atomic propagators also become independent of 
sites: For example, the replacement of Xj — > X leads the 
slave boson propagator Q h °J (r 3 ■, ioj n ) Gb{iu n ) where 
the site-independent propagator is defined as 



1 



(8) 



Let us solve the self-consistent equation ([7]). For 
ksT/U <C 1, Eq. ([7]) can be simplified, and solved an- 
alytically as follows: 



3 + ^9 + 8(e-( c/ -^)/^ T + e-"/fcBT) 
A = fc s Tlog - . (9) 

As discussed below, Eq. ([7|) is numerically solved, and we 
compare the numerical result with Eq. (j9|). The obtained 
A leads the estimation of the slave particle densities. The 
density of each slave particle is given by the MF solution 
A: = 6(eb) for a slave boson, n^ F = /(eh) for a 

holon and n^ F = /(e<i) for a doublon. The temperature 
and chemical potential dependency of them are shown 
in Fig. [TJ In the temperature region shown in Fig. [TJ 
the results analytically given by Eq. © are in precise 
agreement with ones given by the numerically solved A. 
It means that in such a regime we may always employ 
Eq. j9} as a solution of the MF self-consistent equa- 
tion ([7|). On the other hand, the densities in the atomic 
limit can be exactly calculated as shown in Appendix 
Here we take the exact densities to be n CT , nh and nd, and 
their temperature and chemical potential dependency is 
also shown in Fig.[T] The exact result allows us to discuss 
the temperature and chemical potential regime justifying 
the slave particle approach. From the comparison made 
in Fig. [TJ the slave particle technique is found to agree 
with the exact result in the temperature and chemical po- 
tential region where n„ = 1/2 and n^ F = = 0. In 
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FIG. 1. (cClor online) The slave-particle densities as a function of temperature for different chemical potentials: (a) the slave 
boson, (b) holon, and (c) doublon density. The chemical potentials are chosen as follows: fj,/U = 0.5, 0.2, 0.1, and from the 
top to bottom. The MF results n^ F , n^ F , and n^ F given by the slave-particle technique and the exact results in the atomic 
limit, n CT , n^, and n^, are compared. nj^ F , nJf F , and n^ F are calculated numerically and analytically. The solid, dashed, and 
dotted lines, respectively, denote the exact result by the atomic limit calculation, the slave particle approach ()A3[) . 



other words the slave particle technique is expected to be 
reasonable when the system is near a Mott insulator, and 
such temperature as a benchmark is below fc^T/f/ k, 0.1. 



B. The finite hopping 

We now consider a finite but small hopping by taking 
the infinite series of diagrams produced by the scattering 
Hk among the slave particles, based on the non-crossing 
(NC) approximation jl3lll4|. This approximation can be 
also regarded as a certain type of high temperature series 
expansion (HTSE) [2(|, but in our formalism the Wick 
theorem is still applicable due to the quadratic Hamilto- 
nian (|S|), and a particular infinite series of kinetic energy 
perturbation can be taken. Based on the NC approxima- 
tion the self-energy diagrams constructed by full propa- 
gators, shown in Fig. [2j are considered. 

Since one is in a temperature regime higher than the 
anti-ferromagnetic exchange ~ 4 J 2 /U, we apply the spin- 
incoherent assumption to the slave boson propagator. 
Namely, the slave boson propagators in the diagrams of 
Fig. [2] are replaced by the atomic one flE), 



(10) 



The doublon and holon propagators are left full ones: 
namely, by using the Dyson equation the full propagators 
are written as 

g d/h (k,iis n ) = Q~d/ h (k,iv n ) 



(11) 



Unlike the standard MF theory, the dynamical fluctua- 
tion of the slave bosons is retained here, which is neces- 
sary to describe the doublon-holon excitation. The NC 
diagram Fig. [2] (b) couples the self-consistent equations 
of the doublon and holon self-energy. However, the con- 
tribution is negligibly small because it is a far off-shell 
diagram in this case such as Mott state. Consequently 
the self-consistent equations of the self-energy Fig. [2] (a) 
are decoupled and one can obtain in momentum space, 
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v 



(12) 



with W = i/8z6(eb)[&(eb) + ^]J 2 corresponding to a half 
band width for the holon and doublon as we will see be- 



( a ) 



(b) 
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FIG. 2. The NC diagrams giving the doublon self-energy. The 
solid, double dotted, and dashed lines, respectively, denotes 
the full propagators of the doublon, holon, and slave bosons. 
The left diagram (a) describes the scattering between a dou- 
blon and slave boson. The right diagram (b) represents the 
higher energy scattering to a holon than the left (a). Thus, as 
long as charge excitations of energy ~ U are taken for large 
U, the diagram (b) would be irrelevant. The holon self-energy 
is also given by the same type of diagrams. 
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FIG. 3. (Color online) The doublon-holon band width 
W/y/zJ as a function of k B T for fj,/U = 0.5. The chemical 
potential dependence of the band width for ksT/U = 0.01 
(solid), 0.05 (dashed), 0.1 (dashed-dotted) , and 0.5 (dotted) 
is also shown in the inset. 



IV. DOUBLON PRODUCTION RATE 

We calculate the DPR induced by the amplitude mod- 
ulation of an optical lattice potential using the above for- 
malism. The amplitude modulation of the optical lattice 
potential, V(t) = Vo + SV cos(ujt), modifies both J and 
U as J -)• J[l + SJcos(ujt)} and U -> U[l + SU cos(wt)], 
where 5 J and SU are dimensionless modulation pertur- 
bation parameters and given as a function of 5V. How- 
ever, it is possible to map the two parameter modu- 
lations problem to single parameter one of either J or 
U [1, [l(J, [n|, HH . Namely, the modulation perturbation 
to be discussed here can be written as 



H mod (t) = SF cos(ujt)H K 



(16) 



low. z is a coordination number, and N is the total site 
number of the system. Note that due to the momentum 
dependence of the right-hand side the self-energy should 
be given as a local quantity: T,^(k,iu n ) — H^(iv n ). 
Consequently the propagator also turns out to be lo- 
cal: Q ( ] l (k,ii' n ) = Gd(iv n ). Thus the self-consistent equa- 
tion (IT21) is easily solved through the Dyson equation (fTTj) 
as follows: 



^d(^n) = 



- e d - i\/W 2 - {iv n - e d ) 2 



(13) 



where SF = 5 J — SU. [2l| Within the second-order per- 
turbation in terms of SF, the DPR defined as the time 
averaged growth rate of atoms forming doublons is given 
as B 



(17) 



where XkH = -i J™dt e iwt ([H K (t), H K (Q)]}. As shown 
in Appendix [B] we represent the correlation function 
Xk(t) = — (T t Hk(t)Hk{0)) without vertex corrections 



Through the analytic continuation the doublon spec- 
tral function, which is equivalent to the density of state 
(DOS) in this case, is also obtained. In this approxima- 
tion, a semi-circle type DOS is formed: 



v ' wy v w 



(14) 



One can also obtain the self-energy and spectral func- 
tion of the holon in the same way. The forms are the 
same as what is obtained by the replacement e d — > in 
Eqs. (fT5|) and (TT4")) . The holon spectral function is found 
to reasonably reproduce the result of Brinkman and Rice 
on single hole dynamics in a Mott insulator The 
chemical potential dependency of the band width W is 
shown in Fig. [3l 

Using the representation ^fy, via the Matsubara 
Green's function of the original fermion, the DOS (spec- 
tral function) is represented as 



b{e h ) + f{e h -uo) 



A(eb - w) 



b(e-b) + /(eb+w) 



^.d(eb 



(15) 



As expected, the doublon and holon spectral functions, 
Ad and Ah, give the upper and lower Hubbard band, 
respectively. The spectral function as a function of J/U 
and the chemical potential is shown in FigHJ 
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FIG. 4. (Color online) The spectral function of the original 
fermion as a function of k B T ' jU for fi/U = 0.5 and J/U = 0.02 
(top panel), n/U for k B T/U = 0.1, and J/U = 0.02 (middle 
panel) and J/U for k B T/U = 0.1 and fj,/U = 0.5 (bottom 
panel) . 
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FIG. 5. The diagrams contributing to the considered correlation function (|18p : The diagrams (a), (b), (c) and (d) correspond 
to the first, second, third, and fourth terms. The solid, double dotted, and dashed lines denotes the doublon, holon, and slave 
boson propagators, respectively. 



in the slave-particle description 22} , 

Xk(t) = -2J 2 X^Tj&M&C-r) + T%{r)Q*{T)Q h {-T) 

- [G b (-r)} 2 g h (r)g d (r) - [g h (T)} 2 g h (-r)g d (-r)}, 

(18) 

where r|(r) = (T r xJ(T)af,-(T)ast(0)xi(0)> (x = d or fc) 
is a two-particle Green's function of a doublon (X = d) 
and holon (X = h). The diagrams corresponding to the 
terms in Eq. (fT5|) are illustrated in Fig. [5] Without the 
vertex correction, the two-particle propagators are con- 
tracted to a single particle propagators by the Wick ex- 
pansion: r*(r) = — Gx(t)Qx(— t). Through the Fourier 
transform of Xk(t) and analytic continuation, one can 
straightforwardly obtain the real-time kinetic energy cor- 
relation function in frequency domain. As a result, the 
imaginary part of the correlation function is given as 



-NW 2 /8 



dp 
2^ 



f(u - w) - f{v)\ [Ah(v)A h (v - w) 



S-[6(2eh) 

27T 



+ /(*/)] {[/(2e b -*/)-/(2e b -^ + a;)] 
x A&(v)Ah{2th + w - i/) - (w -> -u>)}. 



(19) 
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FIG. 6. (Color online) The DPR spectrum as a function of 
modulation frequency. The solid line and points denote the 
theoretical and experimental results, respectively. The tem- 
perature necessary to draw the theoretical curve is determined 
in the inset; the temperature dependence of the theoretically 
given DPR at u> = U/h and the temperature is determined 
from the crossing point to the experimental data (dotted line), 
which is k B T/U » 0.052. 



our framework, temperature is treated as a free parame- 
ter so that we determine the temperature by a fit of the 
DPR spectrum intensity at uj = U/H, which is obtained 
in the experiment. 



V. COMPARISON WITH EXPERIMENTS 

Let us compare our result (|19p with the experimental 
data. We employ a set of the parameters evaluated in 
the 40 K atom experiment @: the hopping J/h = 2n x 85 
(Hz) and the interaction U/h — 2tt x 5400 (Hz). In terms 
of the optical lattice potential, the depth, modulation 
rate and lattice constant are, respectively, taken to be 
Vb = 10.Er where £r is a recoil energy, SV/Vq = 0.1 and 
a — 532 (nm). The lattice modulation is translated into 
SF R3 —0.32 in hopping modulation. The LDA is used to 
take into account the effect of the harmonic trap poten- 
tial Vtrap(f) whose frequency is (lu x /2tt, uj y /2ir, uj z /2it) = 
(56,61,139) (Hz). In the LDA, we replace the chemi- 
cal potential of the homogeneous case by the local one, 
fj,(r) = fio — Vtrapif) where is self-consistently deter- 
mined to give the total trapped atom number 8 x 10 4 . In 



The temperature dependence of the DPR spectrum at 
lj = U/h is shown in the inset in Fig. |6l and U-qT/U ~ 
0.052 in this system is determined. [23| The determined 
temperature is in the region to justify the slave-particle 
approach, which is discussed in Sec. lIIIl and the theory is 
thus expected to work well. Furthermore, using the ob- 
tained temperature, we plot the DPR spectrum in Fig. [5] 
The agreement is remarkably good. In addition to giv- 
ing access to the line shape it means that via our theory 
one can use the shaking method as a good thermometer 
in low temperature regime |24| . since the curve giving 
the amplitude versus temperature (the inset in Fig. [5]) 
is reasonably smooth and steep. To check this point we 
compare in our case the temperature determined by the 
fitting of the shaking curve with other estimates from en- 
tropy calculations [25[ and find that the two results are 
perfectly consistent. 
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VI. CONCLUSION AND SUMMARY 

In this paper, we have described the charge excita- 
tion of strongly correlated fermionic systems in the spin- 
incoherent paramagnetic regime by a slave-particle rep- 
resentation and diagrammatic approach from the atomic 
limit. This method allows us to take the finite tempera- 
ture and trapping into account. Based on the spectrum 
functions of the doublon and holon, the analytical form 
of the DPR spectrum as a second-order response of the 
optical lattice modulation has been given, and extended 
to the homogeneous system of the trapped atom cloud by 
using the local density approximation. In addition, sub- 
stituting the parameters evaluated in the experiment Q , 
a comparison with the experiment has been made. Al- 
though temperature has not been directly measured in 
experiments, it has been determined as an optimization 
parameter which is controlled to fit the experimental 
data. The result has been in agreement with the exper- 
imental data, which shows that one can use the lattice 
modulation spectroscopy as a thermometer. 

Our method has potential extension such as SU(iV) 
higher symmetric atom systems realized in alkaline- 
earth-metal atom experiments. (26l - [29j 
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Appendix A: The slave-particle density in the 
atomic limit 

In order to clarify the parameter regime to justify the 
slave-particle approach, we make a comparison with the 
exact calculation in the atomic limit. Then we focus on 
the single site occupations as a function of temperature 
and chemical potential, which is identical to the densities 
of the slave boson, holon, and doublon. 

Let us perform the atomic limit calculation. We con- 
sider the atomic Hamiltonian (|lb[) . Then the partition 
function is calculated as 



■^atom 
*H/K B T 



l + 2x 



x 2 y 



(Al) 



where x = e^ /KBl and y — e~ u / k]BT . To calculate 
the densities of the slave boson, holon, and doublon, 
we introduce the projection operators. In the original 



fermion picture, the projection operators of the slave bo- 
son, holon and doublon state are, respectively, written as 
Pa = n a - n-f-n^, ph = (n-|- - l)(nj. - 1), and p d = n\n\,. 
They obviously obeys the constraint ~^2 a Pa+P\i+Pd = 1- 
Taking the statistical average of these projection opera- 
tors, we obtain the following exact result of slave-particle 
densities: 

x 

n a = (Pa) = -—: — ; — 5-, (A2a) 



(Ph) 



«d = (Pd) 



1 + 2x + x 2 y ' 
1 

1 + 2x + x 2 y ' 

x 2 y 
1 + 2x + x 2 y 



(A2b) 
(A2c) 



On the other hand, in the case of the slave-particle 
approach combined with the MF solution of the Lagrange 
multiplier A, the slave-particle densities are given as 

„MF 



MF 



MF 



= b(e b ) 



(A3a) 
(A3b) 
(A3c) 



where eb = A, eh = H + A, and ed = U — fi + A, and 
A is a solution of the MF constraint equation ([7]). For 
ksT/U <C 1 the MF Lagrange multiplier A is analytically 
obtained as shown in Eq. The comparison of the 

exact result (| A2|) and the MF approach (| A3|) is shown in 

Fig.ru 



Appendix B: Calculation of the kinetic-energy 
correlation function 



Here we show the calculation of the kinetic energy cor- 
relation function xk(t) in detail, and derive Eq. (|18p 
given in Sec. IIV1 Using the expression of the kinetic en- 
ergy Hamiltonian (TJ|, the correlation function xk(t) is 
straightforwardly written as 



XK 



{t) = -J 2 E {{TrF 3l {r)F m m){T T [h\{r)h 3 {r) 



(i,j),{l,m) 



- d\{T)d 3 (r)][h\{Q)h m {Q) - 4(0)d m (0)}) 

+ {T T A\. (r)A lm (0)) {T T dj (t)<& (0))(T T hi {r)h\ (0)) 

+ (TrAii {r)A\ m (0)}(T T h\ {r)h m (0)> (T T d] (r)d, (0)) }, 

(Bl) 

where the correlation functions between the operators 
Fji and A m i have been supposed to be zero because the 
slave bosons do not condense in this case. The autocor- 
relations for Fji and A m i should be finite, and in the cal- 
culation the spin-incoherent assumption (|10[) is applied 
in the same way as in Sec. Mil As a result, the autocor- 
relations are calculated as follows: 

(T T F Jt ( T )F lm (0)} = 25 jtl S id g h (-T)g h (T), (B2) 

{T T A VJ (T)A\ m {Q)) = -2 foAn. + <Wi,0 & (t)&(t). 

(B3) 
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As seen in Sec. IIIII the doublon and holon propagators We substitute Eqs. (|B2|) - (IB5|) into Eq. (|B1[) . and then 
turn out to be local: Namely, the correlation function is rewritten as 

Xk(t) = -2J 2 £{(?b(-T)C7b(r) [I* (t) + r^-(r)] 

- [a b (-r)] 2 g h (r)(? d (r) - [<? b (T)] 2 C7 h (-T)C7 d (- r )}, 

(B6) 

where the following two particle correlation functions 
have been defined as 

I* -(t) = (T^t^t^MO)), (B7) 
(T T h i (T)hU0)) = iSiiQuM, (B4) , , t 

{T r di(T)4(0))=t^g d (r). (B5) 
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